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Improving the understanding of strongly correlated quantum many body systems such as gases of interacting 
atoms or electrons is one of the most important challenges in modern condensed matter physics, materials re¬ 
search and chemistry. Enormous progress has been made in the past decades in developing both classical and 
quantum approaches to calculate, simulate and experimentally probe the properties of such systems. In this 
work we use a combination of classical and quantum methods to experimentally explore the properties of an 
interacting quantum gas by creating experimental realizations of continuous matrix product states - a class of 
states which has proven extremely powerful as a variational ansatz for numerical simulations. By systematically 
preparing and probing these states using a circuit quantum electrodynamics (cQED) system we experimentally 
determine a good approximation to the ground-state wave function of the Lieb-Liniger Hamiltonian, which 
describes an interacting Bose gas in one dimension. Since the simulated Hamiltonian is encoded in the mea¬ 
surement observable rather than the controlled quantum system, this approach has the potential to apply to exotic 
models involving multicomponent interacting fields. Our findings also hint at the possibility of experimentally 
exploring general properties of matrix product states and entanglement theory. The scheme presented here is 
applicable to a broad range of systems exploiting strong and tunable light-matter interactions. 


Progress in revealing relations between solid state physics 
and quantum information theory has constantly extended the 
range of quantum many body problems which are tractable 
with classical computers. One such successful approach is the 
density matrix renormalization group (DMRG), which was in¬ 
troduced by White in 1992 [1] and since then developed into 
a leading method for numerical studies of strongly interacting 
one dimensional lattice systems [2]. Later it was realized that 
the DMRG can be interpreted as a variational optimization 
over matrix product states (MPS) [3, 4]. The class of matrix 
product states [5, 6] naturally incorporates an area law for the 
entanglement entropy [7] and is thus ideally suited to parame¬ 
terize many-body states with finite correlation length [8]. An 
interesting connection between the MPS formalism and open 
quantum systems has recently been discovered [9], which led 
to the suggestion of using the high-level of experimental con¬ 
trol achievable over open cavity QED systems to create con¬ 
tinuous matrix product states [10] as itinerant radiation fields 
for the purpose of quantum simulations [11, 12]. In this letter 
we provide experimental evidence that this concept is indeed 
capable of determining the properties of strongly correlated 
quantum systems and offers promising perspectives, comple¬ 
mentary to existing digital and analog quantum simulation ap¬ 
proaches [13] explored with trapped atoms and ions [14-18]. 


In clear distinction to previous experiments, here, we sim¬ 
ulate a continuous quantum field theory rather than a lat¬ 
tice model. In particular, we study the ground-state of the 
Lieb-Liniger model Tf = f dxH = f dx(]V -b T + W) 
[19], which describes a gas of interacting bosons confined 
in a one-dimensional continuum [20], as schematically de¬ 
picted in Lig. la. Here, N = —fiipl.'ipx is the potential en¬ 
ergy, f = da:'tplda;1px is the kinetic energy of particles, and 
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EIG. 1: Schematic of the interacting Bose gas and the principle of 
the quantum variational algorithm, (a). Bosonic particles with ki¬ 
netic energy T are propagating in one dimension along the x axis. 
Repulsion between particles mediates an interaction energy W. The 
particle density p of the gas is controlled by the chemical potential 
p. (b). We experimentally simulate the ground-state of H hy em¬ 
ploying a variational minimization procedure. A tunable cavity QED 
system is used to generate radiation fields emulating continuous ma¬ 
trix product states |0(A)) in a ID transmission line. The average 
energy E\ = (<^(A)|77|(^(A)) of the simulated Hamiltonian is ex¬ 
perimentally determined from measured correlation functions. Ex¬ 
ternal control fields A are used as variational parameters. (c),(d), Ex¬ 
amples of first-order and second-order correlation 

functions measured (dots) in superconducting circuits and simulated 
(solid lines) using a master equation approach for the indicated drive 
rates fl, effective anharmonicity cxj^-K = 5.2 MHz, and cavity decay 
rate K/27r = 2.2 MHz. 
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FIG. 2: Measured correlations in the variational space spanned by the anharmonicity a and the drive rate G. (a) Average photon flux (0) 
proportional to the potential energy {N), (b) J AujLiP 'proportional to the kinetic energy (T), and (c) second-order correlator 
proportional to the interaction energy (VF). The smallest drive rate is Gq/Stt = 0.37MHz. 


is the interaction energy expressed in second 
quantization by the field operator ipx- The Lieb-Liniger model 
has only one intensive parameter v = v/p, where p = (ipl.'ipx) 
is the average particle density and v the interaction strength. 
The ground-state energy of this model can be calculated ana¬ 
lytically using the Bethe ansatz [19]. The calculation of two 
point correlation functions requires the use of numerical meth¬ 
ods such as quantum Monte Carlo or DMRG [10]. The fact 
that this model is well understood makes it an ideal test case 
to benchmark the as yet unexplored quantum variational algo¬ 
rithm recently proposed by Barrett et al. [12]. 

In the experiments presented here, we prepare continuous 
matrix product states |</>(A)) [10, 12] as microwave radiation 
helds propagating along a one-dimensional transmission line, 
see Fig. lb. The radiation helds are generated by an ancillary 
quantum system - in our case a tunable circuit QED system 
[21] - which is coupled with rate k to the transmission line. 
Notably, any radiation held generated in this way is described 
by a continuous matrix product state with a bond dimension 
depending on the number of participating ancillary energy lev¬ 
els [9]. We vary the quantum state |(/)(A)) by tuning a set of 
two external variational parameters A = (a, 12). Here, fl is 
the drive rate of a coherent held applied resonantly to the up¬ 
per eigenmode of the coupled system and a is the effective 
anharmonicity of the driven mode. Tunability of the effective 
anharmonicity a is experimentally achieved by employing a 
qubit of which both the frequency and its coupling to the cav¬ 
ity are adjustable in-situ during the experiment [22], see Ap¬ 
pendix A 5 for details. 

The variational ground-state of the Lieb-Liniger Hamilto¬ 
nian H is found by measuring the expectation value E\ = 
(^(A)|iT|0(A)) = {N) + (T) + {W) and minimizing Ex 
with respect to different variational states |^(A)) created in 
our experiment. The simulated Hamiltonian is thus solely 
determined by the measurement observable. Any model of 
which the corresponding expectation values can be measured, 
is therefore accessible with this approach. Given the pho¬ 
tonic realization of |(/>(A)), the measurement of Ex trans¬ 
lates into the measurement of photon correlation functions. 


Spatial cotTelations in the field are mapped onto time 
correlations in the cavity output field aout(f) by identifying 
V'x = flout (f = x/s)/x/s, where the scale parameter s = xjt 
acts as an additional variational parameter [12]. Entanglement 
in the matrix product states thus corresponds to entanglement 
between photons emitted from the cavity at different times. 
According to this correspondence. Ex for the Lieb-Liniger 
Hamiltonian is given by the hrst- and second-order correla¬ 
tion functions = (aQ^j(T)aout(O)) and = 

(a],ut(0)dLt('^)«out('r)aout(0)) [23]. More specifically, the 
average kinetic energy (T) = s~^ f da;w^G*'^^(a;) is cal¬ 
culated from the Lourier transform of the hrst-order corre¬ 
lation function G^^^(w), the interaction energy is {W) = 
s“^uG*'^^(0) and the potential energy is given by the average 
photon flux (TV) = — s“^/rG^^)(0). 

The presented variational approach thus crucially relies on 
the ability to generate and probe a wide range of different 
(quantum) radiation fields with high efficiency. Lor fast and 
reliable correlation measurements [24] we have developed a 
quantum-limited amplifier which allows for phase-preserving 
amplification at large bandwidth and high dynamic range [25]. 
The examples of measured correlation functions shown in 
Lig. Ic-d illustrate their dependence on the drive rate H at 
constant a. While = 0) equals the total average pho¬ 
ton flux (ao^^aout), the limit —>■ cxd) is proportional 

to the square of the coherence of the held | (flout) P- There¬ 
fore, G^^^ increases with drive rate due to the enhanced pho¬ 
ton production rate. The normalized second-order correlation 
functions = G^^)(t)/(G(^)( 0))^ show anti-bunched 

behavior (0) < 1) for weak drive and Rabi type oscilla¬ 
tions when the drive rate H becomes larger than the decay rate 
[24]. Both the measured hrst-order and second-order correla¬ 
tion functions are in agreement with the results obtained from 
master equation simulations (black solid lines). 

We have measured ^ 10^ such correlation functions over 
a wide range of variational parameters f2 and a by acquiring 
data for one week. Without the employed parametric ampli¬ 
fier the time for measuring this set of data would have been 
on the order of years because of the exponential scaling be- 
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tween statistical error and correlation order [26]. Based on 
this collection of time-resolved correlation functions we have 
evaluated the three relevant terms J (oj) 

and that enter the calculation of E\, see Fig. 2. As 

expected, the average photon flux (Fig. 2a) increases 

with drive rate (1 and is suppressed for increasing anharmonic- 
ity a. The kinetic energy term J duJui'^G^^\uj) in panel b is 
determined by the power spectral density Only the 

spectral weight of photons generated at finite detuning from 
the drive frequency (|w| > 0) contributes to the integral. In 
the Bose gas picture such photons correspond to particles in 
finite momentum states and therefore carrying kinetic energy. 
The rate of scattering events from drive photons into photons 
with finite uj increases with drive strength and has a non-trivial 
dependence on a. Finally, the second-order correlator (0) 
in Fig. 2c clearly reveals the crossover from antibunched ra¬ 
diation —>^ 0) for large a and small 17 to coherent 

radiation (0) —?► 1) when lowering the anharmonicity. 

Based on the three measured quantities shown in Fig. 2 and 
a chosen interaction parameter v we evaluate the energy Ex 
for all prepared states |0(A)). We identify a local minimum 
in variational space (blue region), which corresponds to the 
variational ground-state of the Lieb-Liniger model, see Fig. 3. 
When changing the interaction parameter v we find the energy 
Ex to be minimized by a different set of parameters (a, 17) in 



FIG. 3: Measured energy landscape for the Lieb-Liniger model. Ex 
calculated from the measurement data shown in Fig. 2 relative to its 
minimum i?min as a function of a and G. Interaction strength v in¬ 
creases from bottom to top as indicated. The color scale is adjusted in 
each panel such that the maximal value max appears red. The min¬ 
ima from bottom to top are located at (G/Oq)^ G {13,10, 7, 3} dB 
and a/2'K G {1.56, 2.25, 3.43, 5.04} MHz as indicated by red dots. 


variational space. While for large values of v the minimum 
appears in the anti-bunched region (lower right corner in top 
panel), the minimum moves to the region where the radiation 
is mostly classically coherent when weakening the interaction 
strength (upper left corner in bottom panel). The maximum 
and minimum value of v which can be explored in this way 
is practically limited by the range of variational parameters 
a and 17 for which correlation functions have been acquired 
experimentally. 

After identifying the variational ground-state for each in¬ 
teraction strength v as the respective minimum in the energy 
landscape, we further investigate its properties. We com¬ 
pare the experimental results with the numerically exact re¬ 
sults obtained from a variational matrix product state algo¬ 
rithm executed on a classical computer with bond dimension 
17 = 14 [10]. We followed the usual convention and rescaled 
all quantities so that they correspond to a particle density of 
p — 1 [10]. The Lieb-Liniger ground-state energy density 
Ell, which is the sum of kinetic energy and interaction en¬ 
ergy, increases with interaction strength and ideally converges 
towards the Tonks-Giradeau limit lim Ell = tt^/S indicated 

•D—>-oo 

as a dashed line in Fig. 8a. Given the small number of varia¬ 
tional parameters the experimental data (blue dots) reproduces 
the characteristic dependence of Ell on of the exact solu¬ 
tion (red solid line) quite well. 

Importantly, having physical access to the ground-state 
wave functions |^(A)) we can also probe quantities beyond 
the ground-state energy, such as two-point correlation func¬ 
tions. First-order correlation functions {iplipo) are obtained 
from by converting time into spatial coordinates 

(Fig. 8b). As expected, we observe a decrease in correlation 
length with increasing interaction strength v. Due to the ab¬ 
sence of spontaneous symmetry breaking in one dimension 
[27], the exact ground state of the Lieb-Liniger model does 
not exhibit Bose-Einstein condensation. The observed finite 
limit lim {ipl,ipQ) is a characteristic feature of matrix product 

x—¥oo 

states which do not support the C/(l) symmetry of the model 
for finite bond dimensions. 

The nontrivial nature of the ground-state in the presence 
of interactions also becomes manifest in the particle-particle 
correlator shown in panel c. With increasing 

V particles are more likely to repel each other leading to anti¬ 
bunching. Our experiments clearly resolve this crossover from 
a weakly into a strongly interacting Bose gas by accessing 
variational wave functions for interaction parameters v over 
two orders of magnitude. While this general behavior is qual¬ 
itatively well reproduced, an accurate quantitative agreement 
with the numerical results would require a larger number of 
independent variational parameters in the experimental real¬ 
ization. This becomes apparent when comparing the experi¬ 
mental results with numerical calculations based on continu¬ 
ous matrix product states of different bond dimensions, D = 2 
and D = 14, where the number of variational parameters is 
2D^. Correlation functions simulated with low bond dimen¬ 
sion D = 2 deviate from the exact results {D = 14) similarly 
as the measured ones (Fig. 8d-g). 

In summary, we have experimentally revealed connections 
between open quantum systems, the matrix product state vari¬ 
ational class, and quantum field theories that can be used for 
practical quantum simulations. The presented quantum varia¬ 
tional algorithm is general in the sense that it can be applied 
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FIG. 4: Comparison between experimental simulation and numerical result, (a) Lieb-Liniger ground-state energy density Ell vs. interaction 
parameter v at constant particle density p — 1 on a log/log scale, (b)-(c) Experimentally obtained first-order and particle-particle 

correlation functions {'iplipl.'ipxipo) for the seven indicated interaction strengths v. (d)-(g) Corresponding numerical solutions using continuous 
matrix product states with hond dimensions D = 2 and D = 14. 


to any one-dimensional quantum field theory. Exploring inter¬ 
acting vector field models seems particularly appealing, since 
they are difficult to simulate on classical computers. Experi¬ 
mentally, this could be achieved by coupling tunable quantum 
systems to multiple transmission lines each representing one 
of the vector field components [28]. Higher accuracy in the 
simulation will require more variational parameters and quan¬ 
tum systems with more degrees of freedom, which is achiev¬ 
able with tunable superconducting circuits. In this context, the 
collective dissipation of multiple emitters coupled to the same 
transmission line at finite distance [29] may turn out very use¬ 
ful. Extensions of the presented approach may be envisaged to 
also explore dynamical phenomena and discrete lattice mod¬ 
els using experimentally created matrix product states. 
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Appendix A: Experimental details 

1. Measurement setup, sample fabrication and 
characterization 

The experiments presented in the main text are performed 
with a device consisting of a superconducting circuit, see 
Eig. 5a for details about the experimental setup. The sample 


(Eig. 5b) consists of a A/2 transmission line cavity with a res¬ 
onance frequency a;res/27r « 7.3425 GHz of the fundamental 
mode. The resonator has one output port which dominates 
the total decay rate k/2tt « 2.2 MHz and one weakly cou¬ 
pled input port Kin/K « 0.01 which is used for coherent driv¬ 
ing of the cavity field. The resonator is fabricated using pho¬ 
tolithography and reactive ion etching of a Niobium thin film 
sputtered on a sapphire wafer. We have fabricated a supercon¬ 
ducting qubit (Eig. 5c) close to one end of the resonator. Both 
its transition frequency ujge and coupling strength g to the res¬ 
onator are tunable by varying the magnetic fluxes threading 
the two SQUID loops [22, 30]. Elux control is achieved by a 
combination of a superconducting coil mounted on the back¬ 
side of the sample holder and a local flux line which couples 
predominantly to one of the two SQUIDs. The currents feed¬ 
ing the coil and the flux line are generated by voltage biased 
resistors at room temperature. The qubit is fabricated using 
double-angle evaporation of aluminum on a mask defined by 
electron beam lithography. The qubit decay and dephasing 
times are measured to be Ti « 1.8 /iS and T 2 « 1.2 /is. The 
anharmonicity of the qubit is {oJef — ujge)l2TT « —80 MHz, 
where We/ is the transition frequency from the first excited 
state |e) to the second excited state |/) of the qubit. 

The sample is mounted on the base plate of a dilution re¬ 
frigerator cooled down to a temperature of about 20 mK. The 
qubit and the resonator are coherently driven through atten¬ 
uated charge control lines. The microwave radiation emit¬ 
ted from the cavity is guided through two circulators to a 
Josephson parametric dimer (JPD), which provides quantum- 
limited amplification at large bandwidth and dynamic range 
[25, 31, 32]. A directional coupler is used to apply and in- 
terferometrically cancel the reflected pump field. The ampli- 
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FIG. 5: (a) Schematic of the measurement setup. For details see text, (b) Optical micrograph of the sample. The second qubit gap in the left 
part of the chip is left empty. Enlarged images of the output capacitor (c), of the qubit (d), and of the input capacitor (e) are shown. 


fled signal reflects back from the JPD, passes a bandpass (BP) 
filter, is further amplified by a high electron mobility transis¬ 
tor (HEMT) amplifier, and is down-converted to an interme¬ 
diate frequency (IF) of 25 MHz. After low-pass (LP) filtering 
and IF amplification the down-converted signal is digitized us¬ 
ing analog-to-digital conversion (ADC) and further processed 
with field programmable gate array (FPGA) electronics. 


2. Controlling the qubit frequency and the coupling strength 

We characterize the coupled cavity-qubit system by probing 
the transmission coefficient of the cavity and fitting the data 
to the absolute square of the expression 

t = - -(Al) 

*(Wres — w) -F i(^uig^-u:)+~f/2 + 


which we obtain from input-output theory for the Jaynes- 
Cummings model [33]. Here, w is the probe frequency, 7 is 
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cj/2tT’[GHz] 


FIG. 6: (a) Measurements and fits of the absolute square of the transmission coefficient |f|^ for varying qubit frequency at approximately 
constant coupling strength gj’l'K = 1.75 MHz. Individual data traces are offset from each other by one. (b) Transmission coefficient for 
varying coupling strength g at constant qubit frequency ojge « Ures- 


the qubit decoherence rate and A is a scaling factor. The probe 
power is chosen such that the average number of excitations of 
the coupled resonator qubit system is much smaller than one. 
In this case the qubit may be approximated by an harmonic 
oscillator. We determine the qubit detuning A = uj^es — ^ge 
and its coupling strength g to the cavity by fitting spectro¬ 
scopically obtained transmission data to the above model, see 
Fig. 6. The magnetic fluxes through the qubit SQUID loops 
and with that the qubit parameters g and A are controlled by 
a pair of voltages Vi and V 2 applied to coil bias resistors. For 
small g and A we approximate the relation between (p, A) 
and iVi , V 2 ) by linear equations of the form 

(g\ _ (mil mi2\ (Vi- Ui,o\ 

l^Ay yTO2i W 22 / \V2 —V2,o/ 

We determine the coupling matrix elements rriij and the offset 
voltages Ui 0 and V 2 _o by recording transmission data for pairs 
(Ui, V 2 )k- For each set of data we extract the correspond¬ 
ing parameters (p, A)fc and perform a least-square fit to deter¬ 
mine the model parameters rriij, Ui_o and V 2 p. By inverting 
Eq. (A2) we calculate the voltages Vi and V 2 for a given set of 
desired qubit parameters [g, A). In order to further fine-tune 
the parameters, we have developed an automated calibration 
algorithm which minimizes the deviations from desired tar¬ 
get values by iteratively measuring, fitting and readjusting the 
control parameters Vi and V 2 . With this control procedure we 
are able to independently set the qubit frequency and the in¬ 
teraction strength. We also use this procedure to monitor and 
correct for slow qubit frequency drifts occurring during long 
runs of the experiment. 

We demonstrate individual control of qubit parameters by 
either tuning the qubit frequency for approximately constant 
coupling strength g (Fig. 6a) or by varying g for fixed qubit 
frequency (Fig. 6b). For all sets of data we have turned on the 
JPD amplifier and have divided out its frequency dependent 
gain. For the measurements in Fig. 6b we have kept the qubit 
resonant with the cavity (A = 0) and have varied g. These 
measurements demonstrate the ability to tune the system from 


the fast cavity (k ^ p ^ 7) into the strong coupling regime 

(5 » [22]. 


3. Josephson parametric dimer amplifier 


To measure higher order correlation functions efficiently 
while retaining a high level of linearity we employ a Joseph¬ 
son parametric dimer (JPD) amplifier. For details about the 
operational principles of the JPD we refer the reader to [25]. 
The gain G measured vs. signal frequency / is approximately 
described by a Lorentzian function (Fig. 7). The amplifier 
bandwidth at full width half maximum is approximately 35 
MHz. In order to increase dynamic range, we have chosen a 
moderate maximum gain of 17 dB. A measurement of the gain 
as a function of signal power Pin results in a 1 dB compres¬ 
sion point at about -107 dBm which corresponds to a photon 
flux of 3000 /rs“^, see Fig. 7b. The largest photon flux which 
is generated in the described experiments is less than 50 
The JPD amplifier is thus far away from its compression point 
for all measured correlation functions. The improvement in 
detection efficiency becomes apparent when measuring the 
noise power spectral density Ss in units of photons per Hz 
per second when the JPD amplifier is turned ON and when it 
is turned OFF. The effective noise level, which is referenced 
back to the input of the JPD amplifier is decreased by more 
than an order of magnitude when it is turned on. The scaling 
of Ss is based on a comparison between the frequency depen¬ 
dent gain Gs and the JPD amplifier noise [34]. The deviation 
from the quantum limit is due to the additional noise from the 
following HEMT amplifier which is comparable to the noise 
at the output of the JPD. The equivalent detection efficiency 
of the amplification chain is Pamp = ^/Ss « 50%. 
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FIG. 7: (a) Measured gain of the JPD amplifier (blue dots) and a Lorentzian fit (solid red line), (b) Measured JPD gain as a function of signal 
power and equivalent photon flux, (c) Measured noise power spectral density Ss referenced back to the input of the JPD amplifier in units of 
photons per Hz per second when the JPD is turned on (blue) and when it is turned off (red). The detuning S is relative to the frequency of 
a coherent test tone applied to the JPD input at frequency 7.35 GHz. The dashed line indicates the quantum noise limit for phase-preserving 
amplification and detection. 


4. Calibration of drive rate and output power 

We have calibrated the total gain of the detection chain in¬ 
cluding all cable losses in order to reference the measured 
photon flux back to the output of the cavity. Calibrating the 
total gain of the detection chain is equivalent to calibrating its 
detection efficiency ptot- The detection efficiency is typically 
limited by losses between the cavity and the first amplifier 
iVioss) and by noise added during the amplification process 
(Pamp), see Fig. 8. We compare the nonlinear response of the 



coupled cavity-qubit system with master equation simulations 
to perform this calibration. We bias the qubit with parameters 
(A, g) /27r = (3, 5.7) MHz and apply a drive held to the input 
port of the cavity at rate H and resonant with the frequency 
10 + = ojge + 1^12 -F + is? /A = 27r x 7.35 GHz of the 
upper Jaynes-Cummings doublet state. For these settings we 
measure the coherent photon flux k|(o)P and the total photon 
flux k| (a^a) p emitted from the cavity for different drive rates 
H. We fit these data sets to the results obtained from master 
equation simulations leaving the the total gain factor of the 
measurement chain and the absolute drive rate incident to the 
sample as free parameters (Fig. 8b). The equivalent detection 
efficiency resulting from this fit is equal to the inverse of the 
scaled noise level and found to be ptot = t7ioss??amp = 0.27. 
Together with the estimate for the efficiency of the amplifi¬ 
cation chain Pamp stated in the previous section, we extract a 
radiation loss of 1 — pioss = 0.46 between the cavity and the 
JPD, which is reasonable given the components and cables 
connecting the two stages. 



FIG. 8: (a) The detection efficiency of the cavity output field is typi¬ 
cally limited by radiation loss, schematically represented as a beam¬ 
splitter with finite transmittivity pioss, and by noise Ss = gHAp added 
in the amplification chain. The total detection efficiency is given 
by the product ptot = ?7iossPamp. (b) Measurement (points) and 
fit (solid line) of the cavity photon number and absolute square of 
the coherent amplitude for the coupled cavity-qubit system driven 
through the cavity input port at rate Q. The smallest drive rate is 
Q.o/2-k = 0.37MHz. 


5. Drive scheme and variational parameters 

In the original proposal for simulating the Lieb-Liniger 
model with cavity QED it has been suggested to keep the 
qubit resonant with the cavity and use the qubit drive power 
Slq and the coupling strength g as two variational parameters. 
We have experimentally realized this scheme and found that in 
the limit of small coupling strengths g the total emission rate 
becomes extremely small which in turn limits the signal to 
noise ratio. This is because the effective emission bandwidth 
scales like g'^ /k when g < k and thus decreases quadratically 
with g. We have therefore developed an alternative scheme 
for which the photon emission rate remains proportional to k 
even in the limit of small g. We have therefore made use of 
the ability to tune both the qubit frequency and the coupling 
strength. Rather than keeping the qubit at fixed frequency we 
adjust for each value of g the detuning A such that lo+ re¬ 
mains at constant frequency resonant with the drive frequency 
00 + = LOd = 2Tr X 7.35 GHz. The spectroscopy data for the 
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FIG. 9: (a) Spectroscopy data for all values of g used in the quantum simulation experiment. The blue data points on gray background are 
obtained using resonator transmission measurements. The red data points are obtained using qubit spectroscopy measurements. The relative 
scaling and offsets are adjusted to display the data in the same plot. The inset shows the energy levels |n±) for the Jaynes Cummings model, 
(b) Pairs of (g, A) extracted from the data shown in a using the fitting routine described before, (c) Effective anharmonicity of the upper 
branch of the Jaynes Cummings ladder calculated from the pairs {g, A) shown in b. 


qubit bias points used in the quantum simulation experiment 
are shown in Fig. 9a and demonstrate constant uj+ over the 
entire range of coupling strengths. In order to keep w_|_ con¬ 
stant, we compensate the larger splitting when increasing g by 
tuning the qubit further away from the cavity, see Fig. 9b. 

For this specific tuning scheme we find that the effective 
anharmonicity a of the upper Polariton ladder decreases with 
increasing coupling strength g, as shown in Fig. 9c. To il¬ 
lustrate this effect we show a schematic drawing of the en¬ 
ergy levels |n±) of the Jaynes Cummings model for the res¬ 
onant case (A = 0) and for the case of finite qubit detuning 
A < 0. The inverse proportionality between anharmonicity 
a and coupling strength g illustrates the appearance of anti¬ 
bunched behavior for the small values of g and the observed 
coherent radiation for large g values for this bias scheme, see 
Fig. 2 of the main text. The drive rate of a coherent field 
applied to the cavity input port acts as a second variational 
parameter for the quantum simulation. 


6. Measurement of correlation fnnction and master equation 
simulation 

We employ fast real-time signal processing performed with 
an FPGA at a clock rate of 100 MHz for the measurement 
of time-resolved correlation functions. The cavity output 


field is processed, as described in section lA. After digi¬ 
tization we multiply the sampled voltages with digital sine 
and cosine waves of frequency a;iF/27r = 25 MHz to ob¬ 
tain the quadrature components lit) and Q{t), respectively. 
We then apply an FIR filter with an effective bandwidth of 
F/27r = 10 MHz to the quadratures in time-domain to obtain 
the filtered quadratures I{t) and Q{t), which in the follow¬ 
ing we write as the single complex valued amplitude S{t) = 
I{t) + iQ{t). To measure the first-order correlation function 
in S{t) we take the discrete Fourier transform (T) of M time 
traces Si{t) of length 10.24 /is, multiply with their complex 
conjugate, and average 


r(i)(T) = 


1 

M 


M 

i=l 


Similarly we extract the second-order correlation function by 
calculating the absolute square of S{t) before Fourier trans¬ 
forming 


r(2)(r) = J-1 


^ i=l 


We record each of these quantities with the drive field turned 
on, giving FQ^jr),rQ^(r), and with the drive turned off, 
giving rQpp(r),rQpp(T). To avoid effects due to slow 
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drifts we alternate between all four measurements every 
12.5 /is. As explained in detail in reference [26] and as 
demonstrated in many experiments since then [24, 35, 36], 
we can use these four measurements to extract the corre¬ 
lation functions = K(a^(r)a(0)) and = 

(a^(O)a^ (T)a(T)a(0))/(a^a)^ of the output field of the cavity. 
In these expressions, a (a^) is the annihilation (creation) oper¬ 
ator of the intra-cavity field. For the cases in which the aver¬ 
age photon number is small ((a^a)) <C 1 we find (0) val¬ 
ues which are systematically smaller than the corresponding 
master equation simulation. We attribute this to a weak ther¬ 
mal background radiation during the off measurements which 
we correct for [37]. Good agreement between the measured 
and simulated correlation functions is found when correcting 
for a thermal photon flux of nth ~ 0.03//rs in the detection 
band. 

We compare these measurements with correlation functions 
obtained from master equation simulations. For these simula¬ 
tions we describe the system by the Hamiltonian 

Hsys/h = (Wres “ LOd)a^ a + {ujge - U}d)b^ b 

+ ^{b^)^b^ + g{a^b + ab^) (A3) 

expressed in a frame rotating at the drive frequency uJd/2Tr = 
7.35 GHz. Here, b and b^ are annihilation and creation op¬ 
erators for an excitation of the transmon. In addition to that, 
we account for qubit decay 7, qubit dephasing 7^ and res¬ 
onator emission k with standard Lindblad terms. Simulations 
are run in a Hilbert space including 6 resonator and 3 trans¬ 
mon levels. In order to account for the finite detection band¬ 
width when simulating the second-order correlation function 
we employ the techniques described in [38]. In this approach 
we introduce an ancillary mode c of frequency tOd, which is 
weakly coupled with rate e/27r = 20 kHz to the cavity and 


decays with a rate equal to the detection bandwidth F. The 
second-order correlation function in c is then simulated based 
on the total Liouvillian and taken as an estimate for the filtered 
correlation function of mode a. 


Appendix B: Theoretical aspects and data analysis 

1. Calculation of the Lieh-Liniger energy from correlation 
fnnctions 

The expectation value to be minimized (H) = (T) + (W) + 
(N) is composed of the kinetic energy of the bosons (T), its 
interaction energy (W) and the potential energy {N). Accord¬ 
ing to the correspondence between the field operator and 
the time-dependent radiation field a{t) = y/s'ipx=st, each of 
these expectation values is proportional to a specific measured 
correlation function 

{f) = ^ J 

{W) = 4g(2)(0), 

(TV) = -^G(i)(0) = -pp. (Bl) 

s 

Here, G(^)(u;) = J^[G(^^(t)] is the Fourier transform 

of the first-order correlation function normalized such that 
/ du;G^^)(a;) = G(^^(0). The energy terms in Eq. (Bl) thus 
explicitly depend on the scaling parameter s, which is to be 
treated as an additional variational parameter. We explicitly 
minimize {H) with respect to s by solving ^ {H) = 0, which 
results in 


3 / da;G^^)(w)a;^ 

-t;G(2)(0) + sj (G(2) (0))" t;2 + igG^^) (0) / dwG(i) 


(B2) 


Correspondingly we obtain an explicit expression for E\ = 
{H), which only depends on the measured correlation func¬ 
tions and the model parameters p,, v. We find the variational 
ground state for a given set of model parameters (u,p) by 
minimizing E\ with respect to a and H. 

2. Scaling transformation 

We follow the usual convention and study the Lieb-Liniger 
ground state subject to the constraint that its particle density p 
is equal to one. Using the procedure described in Sec. B 1 we 
find a ground state which generally does not obey this prop¬ 
erty. We therefore apply a scale transformation by adjusting 
the chemical potential p such that p —1. Under the following 
transformation 

(p,w) -)► (p,u) = {p.y^,vy), 


the ground state remains invariant up to a change in the param¬ 
eter s ^ s = s/y, which immediately follows from Eq. (B2). 
Any variational ground state can therefore be transformed into 
another ground state satisfying p = 1, by choosing y appro¬ 
priately. We apply the following procedure to perform this 
scale transformation: 

• Chose parameter v, set p = 1, and find the set of varia¬ 
tional parameters (smin, ^min, ctmin) minimizing E\. 

• Evaluate the particle density p = G* *^^^(0)/Si„in at this 
minimum. 

• Calculate the new chemical potential p = p/p^ and the 
new interaction parameter v = v/p. The new scaling 
parameter becomes Smin = SminP = G^^^(O). 

• The variational parameters (smin, ctmin) specify 
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the variational ground state for the model with interac¬ 
tion strength v and unit particle density. 

Ground states for different interaction parameters are obtained 
by starting the above procedure with a different value for v. 
The Lieb-Liniger energy Ell at interaction strength v (Fig. 4a 
of the main text) is given by 

Ell = {f) + {W) 

^min ^min 

= (G^^'>{0)y^ J dujG‘'^\uj)uj‘^ +vg^^\0). (B3) 

Correlation functions for the Lieb-Liniger model are directly 
obtained from the measured correlation functions by identify¬ 
ing 


{■fpxi’o) = 


= x/Smin) 


G(i)(0) 

and analogously for the second-order correlation function 


(B4) 


where pss is the solution of the matrix equation 

0 = -i[K,p\+RpR^-yn^R^p), (BIO) 

and K = iQ + ^R)R. In order to find the variational mini¬ 
mum of 

E[v,p-Q,R\ = (4/[g,i?]|ff(u,/r)|4/[Q,i?]) 

= tr (I [Q, i?]t [Q, R] + vR^^R^ - pR^R^ p^ 

(BID 

with respect to Q and R a tangent-plane method using the 
time-dependent variational principle (TDVP) in imaginary 
time was exploited [41]. This method proceeds as follows. 
Firstly, v and p are selected. Then a value D as large as pos¬ 
sible is chosen. An initial guess |'I'[(5oj Ro]) for the ground 
state results from random choice of Qq and Rq. Also a toler¬ 
ance T] and a step size <5 is selected. Set j = 0 and perform 
the following sequence of operations until the desired conver¬ 
gence is reached. 


(DiDoDoV’x) = = xlsmin)- (B5) 

3. Numerically exact solution 


1. Calculate the gradient \7z\'i'[Qj, Rj]), where z is a 2D^ 
vector containing the entries of Qj and Rj in, say, lex¬ 
icographic order. Thus Vi|T'[Qj,i?j]) is a vector of 
2L>^ cMPS states. 


The exact solution [19, 39], via the Bethe ansatz, of the 
Lieb-Liniger model at unit density is only possible for one 
specific value of the interaction parameter, namely u = 2. In 
order to calculate properties of the model at unit density for 
other values of v it is necessary to take recourse to numerical 
methods. We exploited a variational method over translation 
invariant continuous matrix product states (cMPS) [10, 11, 
40] 


|vl/[Q,i?])=tr 


^7^ exp 



Q ^I-\- R^ipldx 


Ifl), 


(B6) 

where V exp denotes the path ordering of the argument from 
left to right for increasing values of x. The operators Q and 
R act on an auxiliary space C^, and jn) is the Fock vacuum. 
The variational parameters specifying the cMPS are precisely 
the two D X D matrices Q and R. The Lieb-Liniger hamilto- 
nian H (in the presence of a chemical potential) is comprised 
of three terms H = T -\- W N, and the expectation values 
of these three terms can be readily computed [40] in terms of 
the variational parameters according to 


(4/[Q,i?]|f|vI/[Q,i?]) =tr([Q,i?]t[Q,i?]p,,), (B7) 

{'i>[Q,R]\WmQ,R]) =vtT , and (B8) 

(T'[Q,i?]|f|T'[Q,i?]) = -ptT{R^Rp ,,), (B9) 


2. Calculate the gradient WzE[Qj, Rj] of the energy ex- 
pec ation value. 

3. Calculate the inverse G~^ of the Gram matrix Gz z' = 

4. Set z,+i = z,- -SZX, G:^^,\/z'E[Q,,R,]. 

5. Set j = j 1, unpack Zj+i into the two matrices Qj+i 
and Rj+i, and repeat step (1) until convergence of the 
energy expectation values reaches the prespecified tol¬ 
erance rj. 

After the above algorithm has terminated the cMPS corre¬ 
sponding to unit density (and at the rescaled interaction pa¬ 
rameter) is obtained via the rescaling procedure described in 
the previous section. This method was used to calculate a 
cMPS representation for the Lieb-Liniger ground state across 
a range of interaction parameters from u = 0 to u = 1000. 
The value 79 = 14 was used throughout as the results so ob¬ 
tained are indistinguishable from the known exact solutions at 
V = 0, 2, oo. For the other values of the interaction param¬ 
eter the accumulated errors were estimated and found to be 
negligible. 
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